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The body-centered cubic Coulomb crystal of ions in the presence of a uniform 
magnetic field is studied using the rigid electron background approximation. The 
phonon mode spectra are calculated for a wide range of magnetic field strengths and 
for several orientations of the field in the crystal. The phonon spectra are used to 
calculate the phonon contribution to the crystal energy, entropy, specific heat, Debye- 
Waller factor of ions, and the rms ion displacements from the lattice nodes for a broad 
range of densities, temperatures, chemical compositions, and magnetic fields. Strong 
magnetic field dramatically alters the properties of quantum crystals. The phonon 
specific heat increases by many orders of magnitude. The ion displacements from 
their equilibrium positions become strongly anisotropic. The results can be relevant 
for dusty plasmas, ion plasmas in Penning traps, and especially for the crust of 
magnetars (neutron stars with superstrong magnetic fields B > 10 14 G). The effect 
of the magnetic field on ion displacements in a strongly magnetized neutron star 
crust can suppress the nuclear reaction rates and make them extremely sensitive to 
the magnetic field direction. 

PACS numbers: 26.60.Gj, 67.80.-s, 52.27.Aj, 52.27.Lw, 52.25.Xz 

I. INTRODUCTION 

The model of a crystal of point-like charges immersed in a uniform neutralizing back- 
ground of opposite charge was conceived by Wigner [l| to describe a possible crystallization 
of electrons. These Wigner crystals of electrons have much in common with Coulomb crys- 
tals of ions with the uniform electron background. The model of the Coulomb crystal is 
widely used in different branches of physics, including theory of plasma oscillations (e.g., 
[2I), solid state physics, and works on dusty plasmas and ion plasmas in Penning traps (e.g., 



Moreover, Coulomb crystals of ions, immersed in an almost uniform electron background, 
are formed in the cores of white dwarfs and envelopes of neutron stars. The properties of such 
crystals are important for the structure and evolution of these astrophysical objects (e.g., % 
0]). In particular, the Coulomb crystal heat capacity [6] controls cooling of old white dwarfs 
and is used to determine their ages (e.g., Crystallization of white dwarfs can influence 
;heir pulsation frequencies. It can thus be studied by powerful methods of asteroseismology 
8j. Microscopic properties of Coulomb crystals determine the efficiency of electron-phonon 
scattering in white dwarfs and neutron stars, and, hence, transport properties of their matter 

n n 

(such as electron thermal and electric conductivities, and shear viscosity, e.g., [9|, 110J) as 
well as the neutrino emission in the electron-ion bremsstrahlung process [U|. Many-body 
ion correlations in dense matter produce screening of ion-ion (nucleus-nucleus) Coulomb 
interaction and affect nuclear reaction rates in the thermonuclear burning regime with strong 
plasma screening and in the pycnonuclear burning regime (when the reacting nuclei penetrate 
through the Coulomb barrier owing to zero-point vibrations in crystalline lattice). The 
description of various nuclear burning regimes and observational manifestations of burning 
in white dwarfs and neutron stars are discussed in Refs. [l2, Q] and references therein. 
The manifestations include type la supernova explosions of massive accreting white dwarfs, 
bursts and superbursts, deep crustal heating of accreted matter in neutron stars. 

Since the late 1990s certain astrophysical applications have been requiring a compre- 
hensive study of Coulomb crystals in strong magnetic fields. The topic has become im- 
portant due to the growing observational evidence that some very intriguing astrophysical 



objects, soft-gamma repeaters (SGRs) and anoma 
class of sources called magnetars (see, e.g., [15|, |l6 



ous X-ray pulsars, belong to the same 



17j | and reference therein). These are 



thought to be isolated, sufficiently warm neutron stars with extremely strong magnetic fields 
B > 10 14 G. For instance, the magnetic field of SGR 1806-20, inferred from measurements 
of its spin-down rate, is B ~ 2 x 10 15 G [ljj. Magnetars are observed in all ranges of the 
electromagnetic spectrum. They show powerful quasipersistent X-ray emission, bursts and 
giant bursts with enormous energy release. During giant bursts, one often observes quasi- 
periodic X-ray oscillation, which a r e interpreted (e.,, Q) a 8 vibration, of neutron stars 
(involving torsion vibrations of crystalline neutron star crust). It is likely that the activity 
of magnetars is powered by their superstrong magnetic fields. Thus the magnetars can be 
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viewed as natural laboratories of dense matter in magnetic fields. In order to build adequate 
models of magnetar envelopes and interpret numerous observations, it is crucial to know the 
properties of magnetized Coulomb crystals. The main goal of this paper is to study in detail 
Coulomb crystals in an external uniform magnetic field. (The results reported here were 
partially presented in [3].) 

The Coulomb crystals in question consist of fully ionized ions with charge Ze and mass 
M arranged in a crystal lattice and immersed into the rigid electron background (in this 
case, rigid means unpolarizable or incompressible, i.e. constant and uniform). The effect of 
the magnetic field B on the ion motion can be characterized by the ratio 

b = uj B /u p , (1) 

where 

ZeB l47rZ 2 e 2 n 

are the ion cyclotron frequency and the ion plasma frequency, respectively; n is the ion 
number density, while c is the speed of light. It is expected that the magnetic field modifies 
the properties of the ion crystal at b > 1 (see, however, Figs. [2] and [8] below) . In a strong 
magnetic field the approximation of rigid electron background is a bigger idealization of the 
real situation in neutron star crust matter, since higher densities are required to achieve full 
ionization and suppress the polarizability of the electron background. The higher densities 
(and Up) imply smaller b. However, the effective ion charge approximation turns out to be 



successful for analyzing partially ionized systems (e.g., |2l|). The quantity b (also equal to 



Va/c, where va = B/^Airp is the Alfven velocity, p being the mass density) is, actually, 
independent of the ion charge. Hence, one can consider large b in Coulomb crystals at not 
too high densities having in mind the effective ion charge approximation. Despite that, 
the effect of the polarizability of the compensating electron background has to be studied 
separately. 

The closely related problem of magnetized Wigner crystals of electrons was studied in the 
early 1980s by Usov, Grebenschikov and Ulinich 3| and by Nagai and Fukuyama 23, 24]. 



Usov et al. [22J obtained the equations for crystal oscillation modes, studied qualitatively the 
oscillation spectrum, and diagonalized the Hamiltonian of the crystal for a proper quantum 
description of the oscillations in terms of phonons. In addition, the authors investigated 
asymptotic temperature and magnetic field dependences of the specific heat, rms electron 
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displacement from the respetive lattice site, and magnetic moment of the crystal. They 
obtained a soft phonon mode with the frequency Q oc k 2 near the center of the Brillouin 
zone. It resulted in a rather unusual low-temperature specific heat behavior (T 3//2 ) instead 
of the standard Debye law (T 3 ). Also, the authors mentioned the dependence of the crystal 
energy on the magnetic field direction as well as the increased stability of the crystal due to 
a suppression of electron vibration amplitudes. In addition, Usov et al. considered the 
dielectric function of the crystal and studied the effect of the electromagnetic field induced 
by the electron motion. Their main results were, however, of semi-qualitative nature, limited 
to various extreme cases, whereas the present paper focuses mostly on quantitative results, 
pertaining to Coulomb crystals of ions, with an eye to astrophysical implications. 

Nagai and Fukuyama [23[ calculated phonon spectra of magnetized body-centered cubic 
and face-centered cubic (bcc and fee) Wigner crystals and compared the energies of these 
crystals at zero temperature as a function of the magnetic field. The crystal energy was 
calculated as a sum of the electrostatic (Madelung) energy, and the zero-point energy of 
crystal vibrations. The effect of the anharmonic and exchange terms was neglected. The 
electrostatic energy is independent of the magnetic field (as long as the field does not alter the 
lattice structure). The vibration energy does depend on the field because the field modifies 
phonon modes. At zero field the energy minimum is realized by the bcc structure, both with 
and without zero-point vibrations (in general, it is not true for polarizable background, e.g., 
25|| ) . The authors showed that for a sufficiently strong magnetic field and at relatively high 
densities [r s < 275, where r s = a e /ao is the standard density parameter, a e = (3/47rn e ) 1 / 3 , 
and ao is the Bohr radius, n e being the electron number density] the full energy is minimized 
by the fee structure. It is worth mentioning that the authors did not consider the dependence 
of the zero-point energy on the magnetic field direction and performed all calculations for a 
fixed direction (along one of the high symmetry axes of the crystals). However, this choice 
of the field direction for the bcc lattice was not optimum, as far as the lattice energy was 
concerned. Moreover, the energy gain obtained by choosing the optimum direction is of the 
same order of magnitude as the difference between the zero-point energies of bcc and fee 
lattices. 

Nagai and Fukuyama [24] investigated an analogous structural transition between bcc 
and hexagonal close-packed (hep) electron Wigner lattices. The hep lattice was found ener- 
getically favorable for a sufficiently strong magnetic field at r s < 10700. This result seems 
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more robust since the energy difference between bcc and hep lattices is several times larger 
than the energy gain obtained by choosing the optimum direction of the magnetic field in 



the bcc lattice (the field direction adopted in 24J was the same as in 231] ). 



In addition, Nagai and Fukuyama [24j described the behavior of all 6 phonon modes 
of the hep lattice in the magnetic field. Also, they analyzed quantitatively the behavior 
of transverse and longitudinal electron displacements from the equilibrium positions for bcc 
and hep lattices (at zero temperature) and concluded (qualitatively) that in strong magnetic 
fields the crystals became significantly more stable in the transverse direction and somewhat 
more stable in the longitudinal direction. 

The present paper is organized as follows. Section [Til discusses equations for Coulomb 
crystal oscillations in the magnetic field. Section II I II focuses on the phonon spectrum of 
such a crystal with a simple lattice (i.e., one ion in the primitive cell). The Hamiltonian 
of the system is diagonalized in Sec. IIV| which allows one to express the ion displacement 
operator via phonon creation and annihilation operators. Section |V] presents numerical 
calculations of thermodynamic functions of the Coulomb bcc crystal for a wide range of 
densities, temperatures and magnetic fields. In Sec. |VTJ these results are applied to the real 
physical system found in magnetar crust. Section fVIII is devoted to an analysis of the rms 
amplitudes of ion vibrations and to the Debye- Waller factor of the magnetized Coulomb 
crystal. Finally, Sec. IVIIII considers the dependence of phonon spectrum moments on the 
magnitude and direction of the magnetic field. 

The results will be parameterized by the quantum parameter 9 = Huj p /T = T p /T (where 
T and T p are temperature and ion plasma temperature, respectively), and by the Coulomb 
coupling parameter T = Z 2 e 2 /(aT). In this (3/47rn) 1//3 is the ion sphere radius. 

The Boltzmann constant is set equal to /c# = 1. 

II. DISPERSION EQUATION IN THE MAGNETIC FIELD 

Let us consider a Coulomb crystal of ions in a uniform magnetic field B. The Lagrangian 
of an ion is 

L= M^ + Ze Av _ z 
2 c 

where v is the ion velocity, A is the vector potential, and Ze <fi is the potential energy of the 
ion (0 being the scalar potential). Choosing the vector potential acting on the i-th ion as 
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Af — [B X Ui] a /2 = e alSl ' (where Ui is the i-th ion displacement), the Lagrangian of 
the ion system can be written as: 

7 N 

L B = L + — 5>,X 

N 

= Lo + ^^s^XnV- (4) 

i=l 

In this case, n is the unit vector along the magnetic field, N is the total number of ions, 
and L is the field-free Lagrangian. 

Introducing the same collective coordinates as at B = 0, 



iVcell 

Hp V MN 



n\'„ = \ i — — E A1 P exp (ik ■ Ri) , (5) 



fc 



^ P = \/ E < ex p H fc • ' ( 6 ) 



iV 

1 

one can rearrange the Lagrangian as 

i = io + yE^/^p- (7) 

In this case, A; is a wavevector in the first Brillouin zone, index 1 enumerates direct lattice 
vectors R\, N cc n is the number of ions in the primitive cell (N cc n = 1 for simple lattices), and 
index p goes from 1 to iV ce n. Summations in Eqs. ([5]) and ([6]) are over all wavevectors in the 
first Brillouin zone and over all direct lattice vectors, respectively. In the thermodynamic 
limit, with A% p being interpreted as generalized coordinates of the system, one obtains the 
following Euler equation: 

K + E V l»<V + "B6 a ^Al p = , (8) 

where T^ pp i{k) is the dynamic matrix of the lattice at B — 0. The Fourier-transform 
^fcp(^) > ^fcp(^) yields the algebraic equation 

- n 2 A« kp + <(*0</ - ^BE a ^n?Al v = , (9) 
;hat can be solved if the ion vibration frequency Q satisfies the dispersion equation below 



22 1 

det \v a v l, - n 2 5 pp ,5 aP - ^ B £ Q7 V} = . (10) 
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Due to the presence of the third term on the left-hand side, Eq. (Q does not represent the 
eigennumber problem for a Hermitian matrix, and the respective polarization vectors are not 
orthogonal. This complicates the reduction of the Hamiltonian to the sum of Hamiltonians 
of independent oscillators. The proper procedure will be discussed in Sec. HV] 



III. PHONON SPECTRUM 



There are 3iV ce n oscillation modes for a given vector k in the first Brillouin zone. The 
requencies of these modes satisfy the generalized Kohn's sum rule f2| s = N ce \\(u 2 + uj 2 b ) 



23]. 



At small k = \k\ the behavior of phonon frequencies is more complex than without the 
field. It depends substantially on the direction of k with respect to B. We restrict ourselves 



to iV cc n = 1. It is possible to obtain the exact asymptotes of Q at small k (cf. Ref. 23]). In 
this limit, the dispersion equation (flQj) can be written as: 



2 



+ (1 + b 2 ) + E B + F = , (11) 
where 

16nE B = 16nE Q - 16nb 2 (k ■ nf 

+ (k ai ) 2 b 2 [(P + 72) + (a + 2 72 )(fc ■ n) 2 

+ (71 - 372) (n^ + n 2 y k 2 y + n 2 z k 2 z )\ , (12) 
k = k/k, while the coefficients Eq and Fq are the same as in the field-free dispersion equation 



26|: 



8ttE = {k ai ) 2 {{(3 + 72) + (71 - 3 72 ) {klk 2 y + k 2 y k\ + klk 2 )] , 
256tt 2 F = (A;a0 4 [(/3 + 7 2) 2 + 2(/3 + 72 )(7i - 3 72 )^ 

+ 3(7i - 3 72 ) 2 ^ z 2 ] • (13) 

The constants a, 3, and 712 characterize the crystal at B = 0. They were defined and 



calculated in Ref. 2g], and were recalculated for bcc and fee Coulomb lattices in Ref. 20]. 
They are reproduced in Table HJ along with the lattice constant a;, for completeness. The 
subscripts x, y, and z refer to the Cartesian coordinate system aligned with the main cube 
of the respective reciprocal lattice. 
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The asymptote of the smallest frequency fii at k — > can be derived by dropping the 
f2 6 -term and choosing the smallest root of the remaining quadratic equation: 



n 2 



E B - jE 2 B -4(l + b 2 )F 



(14) 



oo 2 2(1 + b 2 ) 

At sufficiently small k, Vl\/uj 2 ~ Fq/\Eb\- Then, for the phonons, propagating strictly 
perpendicular to the magnetic field, VL\ oc k, and the phonons are acoustic. If, on the other 
hand, k ■ n ^ 0, then flj oc k 2 at small k, in contrast with the linear dependence at B — 0. 
As the angle between k and B decreases, the quadratic asymptote of Q\ becomes valid for 
wider range of k. Neglecting angular dependences and numerical factors in Eqs. ffl2|) and 



(113p . one can estimate the lowest phonon frequency as VL\/uj p ~ kai/y/2 + b 2 for propagation 
perpendicular to the field, and as Qi/u p ~ k 2 af/b for k ■ n ^ 0. In both cases, at b ^> 1, 
the mode typical energy is inversely proportional to the field strength. 

TABLE I: Dynamic matrix coefficients for bee and fee lattices 





a 


P 


7i 


72 


naf 


bee 


4.1243864 


0.84911538 


-3.4886939 


-1.5915193 


2 


fee 


4.0036483 


1.2376805 


-4.2930041 


-1.71184285 


4 



The asymptote of the biggest frequency ^3 can be found by dropping the last term in 
Eq. (Ilip and choosing the maximum root of the remaining quadratic equation: 

1 



fif 

uj 2 



l + b 2 + v/(l + b 2 ) 2 + AE B 



(15) 



At small k this yields 



to 2 



1 +V + V(l + b2 ) 2 -46 2 (fc ■ to) 



+ 0(k 2 



(16) 



In general, at any k and at b 1, the biggest frequency ^3 ~ cj^- This corresponds to the 
conventional cyclotron ion motion. 

Consider k = 0. Then fi|/Wp = 1 + b 2 for phonons propagating perpendicular to the 
field. From the sum rule, it follows that the intermediate frequency ^2 becomes at k = 0. 
Because Q 2 > ^1, it is clear that Q 2 oc k for k — > 0. If, on the other hand, the phonon 
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FIG. 1: Phonon spectrum of the magnetized bcc lattice for several values of b = LUg/cUp in natural 
(a) and logarithmic (b) scales. Slowly growing mode O/cjp ~ 3 at b = 3 is not shown in panel (a). 
Magnetic field and wavevector k orientations are defined in the text; a is the ion sphere radius. 

propagates along the magnetic field, then at k = one has Vt\/uo 2 = max(l,6 2 ), and, 
hence, Vl\/uj 2 = min(l,6 2 ). Therefore both maximum and intermediate frequency modes 
are optical. As the angle between k and B decreases from 7r/2 to 0, the value of Q2 at 
k — > increases from to min (u p , Ub), whereas the value of Q% at k — > decreases from 
•y/u; 2 + cu B to max (cu p , cu B )- 

Obviously, Eq. ( TTTT) can be solved analytically without neglecting any terms. However, 
at very small k analytical schemes become numerically unstable. The thermodynamic prop- 
erties of Coulomb crystals at low temperatures require integration of functions containing 
the phonon frequencies near the center of the Brillouin zone. Under these conditions and 
especially at high magnetic fields the asymptotes given here become helpful. 

Figure [T] shows the phonon spectrum of the bcc lattice (in natural and logarithmic scales) 
in the direction of k determined by the polar angle 9 = 1.3 tan -1 \[2 and the azimuthal 
angle <fi = 0.97r/4 [angles 9 and <fi are defined with respect to the same Cartesian reference 
frame as the x, y, z-subscripts in Eqs. f lT2|) and ([TBI ]. The magnetic field is parallel to the 
direction from an ion towards one of its nearest neig hbors (9 = tan" 1 >/2, = 7r/4). For 
ka < 0.01, analytic asymptotes discussed above are used; at higher ka, exact calculations 
are performed. One can observe the quadratic dependence of the lowest frequency on k 
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i a frequency ~ u p into an acoustic mode in the vicinity of the Brillouin 
0]). The background polarization effects in a magnetized crystal will 



near the Brillouin zone center. The dependence becomes linear closer to the Brillouin zone 
boundary (at b <C 1). 

The polarization of the electron background in the absence of the magnetic field converts 
the optical mode wit 
zone center (e.g., {25, 

be investigated elsewhere. A similar conversion effect can be expected in this well, 
except that one of the modes must remain optical close to the center of the zone, with 
a frequency Q ~ ujb- This is so, because the coefficient of the f2 4 -term in the dispersion 
equation f TTUT) . that determines the sum of the squares of all frequencies, will continue to 
contain u B , whereas in the field-free case at k — > the f2 4 -term tends to zero for polarizable 
background. 

IV. HAMILTON! AN OF THE COULOMB CRYSTAL IN THE MAGNETIC 

FIELD 

In this section the Hamiltonian of the Coulomb crystal in the magnetic field is represented 



as a sum of Hamiltonians of independent oscillators. The derivation follows that of Ref . |22j , 
but some additional details are provided. Consider the Cartesian reference frame with the 
axes directed along the eigenvectors of the matrix V a ^(k) (again iV ce n = 1). Then the 
Lagrangian ([7]) has the form 

k k 

+ f E^sn^v ( ir ) 

k 

In this case, u^a is the phonon frequency at B — 0. Omitting the constant equilibrium 
electrostatic energy term Uq and turning to the Hamiltonian Hb, one has 

= A + yA^ 7 _,, (18) 
Hb = J2 A % V %- Lb 

k 

= + ( 19 ) 

k k 

where 7r£ = V% — UB£ al3 ^ 1 n l3 A' y _ k /2. In the quantum formalism, A% and V k become operators 
with the usual commutation rules A k ,V k , = ih5 a,3 5kk'- However, the operators 7r£ and 
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ir_ k do not commute at a ^ (3: 7r^,7rf, = ifvjjBS l3ia n < 5y - k - Consequently, the creation 
and annihilation operators, defined in the same way as at B = 0, but with 7r£ in place of 
V k , do not satisfy the required commutation relationships. 

In this situation one constructs the creation and annihilation operators as linear combi- 

n 

nations of the operators of generalized coordinates and momenta (e.g., Ref. [22]), 



4 = + 



(20) 



where a k and (3 k are constant coefficients (and summation over A is implied). They can be 



determined from the equation 
one arrives at 



H B , ail 



hn k al (e.g., Ref. [22]). Using Eqs. (jI5j| and (1201. 



fc "fc 



A _A 



(21) 



or 



x 



kX "fc ) 



tt k al = iu B e a ^a%-i(3% 



(22) 
(23) 



so that 



kX 



nl)a£-in h u B e*> a n?'afi 







(24) 



The system ( 1241) can be solved only if the quantities Q k satisfy the dispersion equation 
(flOl) . i.e., if they are the eigenfrequencies of the oscillation modes with the wavevector fc. 
Therefore, there are three possible solutions for the operator a k , corresponding to the three 
generally different frequencies Q ks , s = 1,2,3: a} ks = a ks {T^ k + ^a;| A ^4^ fc /f2 fcs ). Evidently, 
the equations for the coefficients a ks and a^ fcs coincide. Thus from now on index fc will be 
omitted where possible. 

The solutions of the system (124j) are sets of the cofactors to any of the rows of the matrix 



22] 



iVt s u)Bn z —iQ s UBTi l 



Q% iVt s ujBTf 



iVL s UBn y —iVLgUjBn 1 



(25) 
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(these sets are proportional to each other). Now one has to normalize these solutions together 
with the operators a\, so that [a s , aj] =1. Accordingly, 



<2 S , Qi s , 



a x *a x ,ihw B e x 7A n 7 
1 1 



n 

where Eq. fl24|) was employed. Now consider the following sequence of equations, 

= a x , {({l 2 s - col) a x * - i{l s u B E X ^n<a u s *} 

= c£ a* - 1 W + 0.0,) (27) 



( + ) , (26) 



[where Eq. ( I24j) was used twice]. Thus the right-hand side of Eq. ([26"]) vanishes, if f2 s 7^ fi s /. 
This property is analogous to the orthogonality condition for the polarization vectors at 
B = 0. At s = s', Eq. ( 126]) determines the normalization of the coefficients a x : 

..A „,A*t / ^A 



All the other commutators of the operators a s and a\ vanish automatically. 

Owing to the equality \H B , a^l = ^Oa^ and the normalization relationships obtained 
above, the Hamiltonian ([T91 can be rewritten in the canonical form 

H B = J2 m ks Ul s a ks + ~ J . (29) 

Moreover, with the help of Eqs. ([19]) and (120]) . it is possible to prove the relationships: 

Y^atai'*-at*a£' = 0, (30) 

s 

k ST^ W\U)\i . x x > ^A*< A'\ rAA' / Q1 \ 

h2_^—^—{a s a s +a s a s ) = 6 , (31) 

s s 

hJ2^(a x a x '* + a x *a x ') = 5 XX ' . (32) 

s 

Multiplying a^_ ks by a^*, and a ks by a^, summing over s, subtracting the second expression 
from the first one, and also using (I3T)]) and (13"T]) . one obtains 

*4fc = a *s a-ks + & ks . (33) 
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The above formula is instrumental in calculations of such crystal properties as the Debye- 
Waller factor and the rms ion displacement from a lattice site in the magnetic field. These 
results are reported in Sec. IVHI 



V. PHONON THERMODYNAMIC FUNCTIONS OF THE MAGNETIZED 

COULOMB CRYSTAL 



Phonon thermodynamic functions in the magnetic field are calculated u sing the same 
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291 ]) as in the 



general formulas (e.g., [28]) and numerical integration schemes (e.g., 
field-free case. The phonon free energy (with phonon chemical potential \x = and neglecting 
zero-point contribution) reads 



ks 



1 — exp f 






T J 



dk 



In 



exp 



Ml 



fc.s 



T 



(34) 



'bz(2tt) 

where V is the volume, and the integral is over the first Brillouin zone. The phonon thermal 
energy E and heat capacity C are then given by 

e nn ka /r _ 

4 ^E -r^7^ - (36) 



E 



(35) 



C — —T 



o 2 f 



Q T 2 



IhV 



fc5 smh\htt ks /2T) 

In the magnetic field, there are two new parameters, b = ujb/oj p and the field direction. 
A study of frequency moments of the phonon spectrum in Sec. IVIIII will show that the 
dependence of these moments on the field direction is rather weak. This allows one to 
expect that the dependence of thermodynamic functions on the field direction is also weak. 
Thus, in the present section the consideration is restricted to the field direction, which 
corresponds to the minimum zero-point energy of the crystal. For the bcc lattice, it is the 
direction from a lattice site to one of its closest neighbors. 

The calculated phonon heat capacity per one ion is presented in Fig. [2] in logarithmic and 
linear scales [panels (a) and (b), respectively] as a function of O^ 1 = T/T p for several values 
of b. In Fig. [2fa) one can clearly see the change of the low-temperature asymptote from T 3 
to T 3 / 2 due to the appearance of the soft mode oc k 2 for non-zero B. Also, this easily 
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FIG. 2: Phonon specific heat (per ion) for the bcc lattice as a function of temperature for several 
values of the magnetic field in logarithmic (a) and linear (b) scales. The magnetic field is directed 
towards one of the nearest neighbors. 

excited mode (e.g., at b = 100) is responsible for a relatively high specific heat C/N ~ 1 all 
the way up to 9 ~ 10 3 . At these temperatures the field-free specific heat is already down by 
6 orders of magnitude. 

In strong magnetic fields, higher temperatures are required to achieve the classical regime 
as compared to the field-free case. This is so, because the classical regime occurs when many 
phonons are excited in all modes. In a strong magnetic field, there is the high-frequency 
(cyclotron) mode, Q3 ~ ujb (Sec. IIIII) . Hence, the classical regime is realized if bd <C 1, in 
contrast to the conventional criterion 9 <C 1 at B = 0. This is illustrated in Fig. [2](b). For 
instance, for b = 100 the classical value C/N = 3 is reached only at 9 < 0.01. 

The energies of three phonon modes are spaced far away from each other in strong 
magnetic fields. The minimum frequency Qi oc 1/B, the maximum frequency f2 3 w lub, 
and the intermediate frequency ^2 ~ (Sec. IIIII) . This gives rise to the pronounced 
staircase structure of the heat capacity seen in Fig. EJ^b) at b = 10 and 100. 

This discussion is further illustrated by Fig. [3], where one can assess qualitatively the 
behavior of the specific heat in various domains of the T-B plane. At high temperatures, 
the crystal is classic and the specific heat reaches its maximum value of 3. In strong magnetic 
fields (b ^> 1) when the temperature drops below Hlob, the cyclotron phonon mode is frozen 
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FIG. 3: A sketch of the specific heat behavior in various domains of the T-B plane. The point, 
where all straight lines cross, corresponds to b = 9 = 1. 

out and C/N « 2. As the temperature drops further, below T p , the intermediate mode S7 2 
freezes out in large portions of the Brillouin zone, and the specific heat, now mainly due to 
the fully-excited soft mode, approaches 1. By contrast, in a non- magnetized (6 c 1) crystal, 
the two lower modes deviate from being acoustic only in the very vicinity of the Brillouin 
zone center, while the third one is ~ uo p . Hence, when T drops below T v in such a crystal, 
the Debye law C oc T 3 is recovered. Finally, at any b, when the temperature is so low, that 
only the least energetic mode contains any heat and the oc k 2 law is probed, C oc T 3//2 
behavior results. 

In Fig. Hlphonon thermal energy and entropy S = (E — F)/T are plotted. In the case 
of energy, plotted in the linear scale in Fig. IU(a), one observes again the staircase structure. 
Though less pronounced, it occurs for the same reason as for the specific heat. The classical 
harmonic oscillator limit E = 3NT is naturally reproduced. At low temperatures, the 
magnetic field gives rise to the T 2 dependence of the energy instead of the T 4 field-free 
asymptote. For sufficiently low temperatures, the thermal energy is dominated by the zero- 
point energy, determined by the spectrum moment u\ and discussed in Sec. IVIIII The 
numerical calculations of the entropy are shown in Fig. H^b). In this CclS6 ; clS for the specific 
heat, the familiar field-free T 3 asymptote is replaced by T 3//2 in the quantum magnetized 
crystal. In the classical regime, the entropy depends on temperature logarithmically and is 
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FIG. 4: Phonon energy per ion divided by T (a) and phonon entropy per ion (b) for the bcc lattice 
as functions of temperature for several values of the magnetic field. Inset in panel (b) shows S/N 
in the linear scale at b = 100. The magnetic field is directed towards one of the nearest neighbors. 



insensitive to the magnetic field. This asymptote is discussed in some detail in Sec. IVIIII 
When drawn in linear scale, the entropy, as a function of log 10 T, also shows an atypical 
structure at b ^> 1 [inset in Fig. IU(b)]. In this case, the dependence becomes piece- wise 
linear with several slope changes, corresponding to the sequential excitation of the three 
phonon modes. 

Clearly, magnetized crystal thermodynamics, constructed in this Section, cannot be de- 
scribed by the well-known in solid state physics Debye model (e.g., 28|). On the other 
hand, at 6 3> 1, partial thermodynamic quantities due to the cyclotron mode can be 
represented with high accuracy by the simple Einstein model (e.g., |30j). 



VI. APPLICATION TO MAGNETAR CRUST 



As a practical application of the previous section results, consider fully ionized iron plasma 
in neutron star crust (Sec. [T|). Typical values of 6 and b at B — 10 15 and 10 16 G can be 
estimated from Fig. [5j The curve V = 175 shows the melting line of the classical Coulomb 
crystal at B = 0. The line T = Z 2 e 2 /a represents typical electrostatic energy per ion (and 
separates the regions of a free ion gas above the line, and a strongly coupled ion system 
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below the line). The electron degeneracy temperature is shown by dotted lines marked T F0 , 
T F15 , and T F16 for three values of the magnetic field, -8 = 0, 10 15 , and 10 16 G, respectively. 
Finally, i?zi5 and Ezw are the electron ground-state energies of one-electron iron ion for 
B = 10 15 or 10 16 G, respectively. These energies are calculated using rescaling of equivalent 



hydrogen energies E(Z, B) = Z 2 E(Z = 1, B/Z 2 ) (e.g., 31[), and the fitting formula [32] for 
the energy of the s = state in the hydrogen atom. The assumption of full ionization in the 
degenerate plasma can be used if T F > E z (for a given magnetic field). At T F 3> E z , the 
electron gas is nearly incompressible, and the rigid electron background model is very well 
justified. At T F < Ez the plasma cannot be treated as fully ionized, but even in this case one 
can use present results for estimates by employing the effective ion charge approximation. 
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FIG. 5: Temperature-density diagram for fully ionized 56 Fe matter. F = 175 is the classical crystal 
melting line. T F and Ez mark electron degeneracy temperature and electron-ion binding energy 
(see text for details). Subscripts 0, 15, and 16 refer to the magnetic field values B = 0, 10 15 , and 
10 16 G, respectively. Note, that Eza = 1.1 x 10 8 K. 



Figure E] demonstrates the phonon and electron specific heat (per ion) of fully ionized 56 Fe 
plasma as a function of density (a) and temperature (b) at B = and 2x 10 15 G. The electron 
contribution is calculated using standard formulas for strongly degenerate relativistic Fermi 
gas. For B = 2 x 10 15 G under the displayed conditions the plasma electrons fill the lowest 
Landau level only; the next Landau level would be occupied at p > 10 9 g cm -3 , and the 
plasma would become partly ionized at p < 10 s g cm -3 (cf. Fig. ED- The temperature 
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dependence of the electron specific heat is linear, and near coincidence of two electron lines 
in Fig. E(b) is largely accidental. For the ions to crystallize, the temperature must be below 
the melting temperature, which is ~ 10 8 K in the density range considered (cf. Fig. [5]). 



i — \ — 1 1 1 1 




FIG. 6: Phonon and electron specific heat (per ion) of fully ionized 56 Fe matter versus density 
(a) and temperature (b) at B = and 2 x 10 15 G. For the B ^ case, b6 = 23 in panel (a); 
T p = 6.2 x 10 7 K and b = 1.09 in panel (b). 

The phonon contribution at B — 0, recalculated here, reproduces the results of {(J. 
At B = 2 x 10 15 G the results of the present work are plotted. As seen from Fig. [6l 
phonons dominate electrons in the specific heat in a wide range of temperatures and densities. 
The magnetic field provides an extra boost to the phonon specific heat, especially at lower 
temperatures (in the quantum regime) due to the easily excited soft mode. 



VII. ION VIBRATIONS AND THE DEBYE- WALLER FACTOR IN THE 

MAGNETIZED CRYSTAL 



Using Eqs. ([5]) and (133]) . one can derive the expression for the operator of the ion dis- 
placement from its equilibrium lattice position: 

ih 



^2{aL a ks - aL* al hs ) exp (ik ■ R x ) . (37) 
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Then the rms ion displacement = a/ (u x u x )t can be calculated as 

r 2 h 2 

ks 

^{hu p cfccfc){2n h . + l) = ~r{e,b) (38) 



3 at 

ks 

[cf. Eq. (EE}], where n ks = (e m "^ T - l)" 1 is the mean number of phonons in a mode ks. In 
addition, it is possible to find the rms ion displacement u q in the direction along an arbitrary 
unit vector q, 

= WnY, a *s < <f <f (2n fcs + 1) • (39) 

ks 

In the field- free case, u q is isotropic (in the bcc crystal), but in the magnetic field it becomes 
anisotropic. 

Amplitudes of ion oscillations are depicted in Fig. [7J The ion equilibrium position is at the 
origin. The F-axis is parallel to B and to the direction towards one of the nearest neighbors. 
There is symmetry in the plane azimuthal with respect to B (the plane perpendicular to the 
F-axis). Hence, the X-axis indicates an arbitrary direction orthogonal to B. The distance 
between the origin and a point on a chosen curve is equal to the rms displacement u q in 
the given direction. The distance is expressed in units of the spacing between the nearest 
neighbors, a n = ai a/3/2 = (37r 2 ) 1 / 6 a (for bcc lattice). Solid, dash-dotted, long-dashed, 
dotted, and short-dashed curves correspond to b = 100,10,1,0.1, and 0.01, respectively. 
The 6 panels of Fig. [7J are for the 6 values of the quantum parameter 9, from 0.1 to 10 4 . 
The ion displacement at a different T (for given 9 and b) can be found by straightforward 
scaling per Eq. (1381) . 

The quantity u q is related to the Debye- Waller factor W(q): (exp(iq • u))t = exp(— W). 
The thermal averaging yields 



(exp (iq • u))t = exp 



MN 

ks 



J2^ a £<sUks + \) 



Thus, u\ = 2W(q). 

As shown in Fig. [7J the ion displacements at 9 < 1 are insensitive to the magnetic field (all 
5 lines merge) and isotropic. At 9 ~ 10 the displacement along the field remains essentially 
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FIG. 7: The rms ion displacement from the equilibrium position (the origin). Y-axis (u q \) is 
parallel, X-axis (it q t) is perpendicular to B. u q \ and u q t are in units of the distance to the nearest 
neighbor. The solid, dash-dotted, long-dashed, dotted, and short-dashed curves correspond to 
b = 100, 10, 1, 0.1, and 0.01, respectively. The distance between the origin and a point on a chosen 
curve is equal to the rms ion displacement u q in that direction. 

the same as at B = 0, whereas the transverse displacements shrink. The net effect will be an 
increase of the crystal stability and melting temperature. The anisotropy of displacements 
in a quantum, strongly magnetized crystal, at 9 > 100, b S> 1, becomes much larger. The 
ion displacements are str ong ly suppressed in the transverse direction and are mildly reduced 
along the field (see also \2^). 

These effects are important for accurate computations of nuclear reaction rates in a 
strongly magnetized crust of a neutron star (Sec. [I]). The strongest effect is expected to occur 
at sufficiently high densities and not too high temperatures in the pycnonuclear burning 
regime. In that regime (e.g., I2I, Q, Q]) the main contribution to the reaction rate comes 
from the closest neighbors in a crystalline lattice due to their zero-point vibrations. In the 
field-free case, the rate depends exponentially on the squared ratio of the equilibrium distance 
between the closest neighbors and the average amplitude of their displacements. The ratio is 



21 



typically large, and the reaction rates are exponentially small (but increase with the growth 
of stellar matter density). Although the reaction rates in magnetized crystals have not been 
studied yet, one can expect a similar result, involving now the average displacement (1391) 
in the direction of closest neighbors. If so, the rates would become extremely sensitive to 
the orientation of the magnetic field. If the field is not directed towards one of the nearest 
neighbors, the reactions will be significantly slowed down by greater tunnelling lengths. Even 
if the magnetic field is directed towards one of the nearest neighbors then, in addition to 
the reduction of the longitudinal displacement, there will be a geometrical effect. In a bcc 
lattice every ion has 8 nearest neighbors, of which only 2 will be located along the magnetic 
field line. The reactions with the other 6 neighbors will be quenched. 

VIII. PHONON SPECTRUM MOMENTS IN THE MAGNETIC FIELD 

According to the well-known Bohr-van Leeuwen theorem, the classical partition function 
is not affected by the magnetic field. For instance, the classical asymptote of the entropy S 
has the form: 



Therefore, the quantity (In (f2/a;p))ph, where average over all phonon modes in the first 
Brillouin zone is implied, should not depend on the magnetic field. This statement is easy 
to prove directly. By inspecting Eq. fflOl) one concludes that its f2°-term cannot contain B 
at any k. On the other hand, this term is equal to f^f^^i) which proves the point. 

Phonon frequency moments u m = ((Q/u p ) m ) ph = u m (b,n) depend on both the field 
strength and direction n. The moment U-2 diverges, while the moment = 1 + b 2 . Consider 
the moments u\ and U-\. Their dependences on the field strength are shown in Fig. [5(a) 
for the field direction corresponding to the minimum of the zero-point energy. At strong 
magnetic fields both moments behave in the same way, proportional to b, although under 
the effect of different modes. The main contribution to u\ comes from the mode 0,% ~ ujb, 
whereas the major contribution to w_i is due to the mode Qi oc 1/B (near the center of the 
Brillouin zone). 

The dependence of the phonon frequency moments on the magnetic field direction turns 
out to be rather weak. In Fig. [HJ(b) the differences Aui = Ui(b, n 2 ) — U\(b, n x ) and Au^i = 
u_i(b, n 2 ) — w_i(6, nx) are shown, where n x is the magnetic field direction towards one of the 




(40) 
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FIG. 8: (a): The dependence of phonon spectrum moments u\ and U-i on the strength of the 
magnetic field directed towards one of the nearest neighbors, (b): The dependence of phonon 
spectrum moment differences Aui and Au_i on the magnetic field (see explanation in the text). 

nearest neighbors, while ri2 is the direction towards one of the next order nearest neighbors. 

The behavior of Au\ is of special importance because it determines zero-point energy 
gain AE zp = ^Nf\w p Aui resulting from the crystal rotation with respect to the magnetic 
field. In weak fields 10~ 3 < b < 1, the difference Au\ scales approximately as b A . In strong 
fields, Am saturates at a value ~ 10~ 2 , depending on the field direction. The difference 
Au-i scales as b 2 for low magnetic fields, and as b for strong magnetic fields. 



IX. CONCLUSION 



The Coulomb crystal of ions with incompressible charge compensating background of 
electrons in constant uniform magnetic field has been studied. The phonon mode spec- 
trum of the crystal with bcc lattice has been calculated for a wide range of magnetic field 
strengths and orientations (Fig. [1]). The phonon spectrum has been used in 3D numerical 
integrations over the first Brillouin zone for a detailed quantitative analysis of the phonon 
contribution to the crystal thermodynamic functions, Debye- Waller factor of ions, and the 
rms ion displacements from the lattice nodes for a broad range of densities, temperatures, 
chemical compositions, and magnetic fields. 
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The characteristic parameter that determines the strength of the magnetic field effects in 
the crystal is the ratio of the ion cyclotron frequency to the ion plasma frequency, b = ujb/oj p . 
Even a very strong field (b ^> 1) does not alter the partition function (hence, thermodynamic 
functions, melting etc.) of a classical crystal, i.e. a crystal at a temperature significantly 
exceeding the energy of any phonon mode (T 3> Hub for strongly magnetized crystal). 

Strong magnetic field dramatically changes various properties of quantum crystals, espe- 
cially at 9 = fiwp/T ^> 1. Low-temperature phonon specific heat, entropy, thermal energy 
increase by orders of magnitude (e.g., Figs. [2, HD- The thermodynamic functions exhibit pe- 
culiar staircase structures, brought about by the vastly different energy scales of the crystal 
phonon mode. Ion displacements from the equilibrium positions become strongly anisotropic 
(Fig. ED, and so does the Debye- Waller factor. 

These results can be applied directly to real physical systems found in the crust of mag- 
netars (neutron stars with superstrong magnetic field). In particular, the heat capacity of 
the magnetized Coulomb crystal composed of fully ionized 56 Fe has been analyzed (Sec. IVII ). 
It has been shown that in the magnetic field (as in the field-free case) the phonon heat 
capacity dominates the electron one in a wide range of densities and temperatures (Fig. ED- 
In accordance with Sec. |V] magnetic field is responsible for a significant boost of phonon 
heat capacity, which may be important for simulations of magnetar cooling. Since ion dis- 
placements are suppressed, one can expect increased crystal stability in the magnetar crust 
associated with an increase of the melting temperature in the quantum regime. Finally, one 
expects a reduction of the pycnonuclear reaction rates due to increased tunnelling lengths. 
The rates are likely to become very sensitive to the orientation of the magnetic field (Sec. 

EH}. 

The dependence of thermodynamic functions on the magnetic field orientation within 
the crystal turns out to be weak, at least for the bcc lattice, but finite. The energy gain, 
achieved by reorienting the crystal with respect to the magnetic field, can be estimated using 
Fig. El 

The techniques used in the present paper can be easily reformulated for fee lattice, where 
results are expected to be numerically close to those for bcc lattice. 
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